/**
* @file curve.cpp
* August 1, 2008
* @author Andre Pereira (brabo@tecgraf.puc-rio.br)
*/

// includes
//----------
#include "curves/curve.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <vector>
#include <set>
#include <algorithm>
#include "curves/line.h"
#include "curves/polyline.h"

using namespace std;

//
/////////////////////////////////////////////////////////////////////
Curve::Curve()
: m_id( -1 )          // Curve id
, m_type( UNDEF )        // type (sub-class)
, m_subdv( 1 )      // subdivisions
, m_ratio( 1.0 )    // subdivision ratio
, m_attr( NULL )      // attributes
, m_pDrawObj( NULL )
, m_firstDraw( true )
{
}


//
/////////////////////////////////////////////////////////////////////
Curve::~Curve()
{
}


//
/////////////////////////////////////////////////////////////////////
void
Curve::getSubdvision( vector<Point3D>& _subdp )
{
  getSubdvision( m_subdv, m_ratio, _subdp );
}

//
/////////////////////////////////////////////////////////////////////
void
Curve::getTangent( vector<Point3D>& _tang )
{
  getTangent( m_subdv, m_ratio, _tang );
}



/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////
//  private methods
/////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////


// given the number of subdvisions and ratio, return 'pa' vector that contains the 
// subdvision in parametric space (0-1). pa.lenght = n+1
/////////////////////////////////////////////////////////////////////
void
Curve::getParamSubdivion( int _n,
                          double _r,
                          vector<double>& _pa )
{
  int i;
  double a;

  // resize pa
  _pa.resize( _n+1 );

  if (fabs(_r-1.0) < 1.0e-4) // equal Lenght 
  {
    for (i = 0; i <= _n; ++i)
      _pa[i] = (double)i / (double)_n;
    _pa[_n] = 1.0; // just to avoid aproximation error
    return;
  }

  a = 2.0 / (_n * (_n-1) * (1.0+2.0 / (1.0/_r-1.0) ));
  _pa[0] = 0.0;
  _pa[1] = a * (_n-1) / (1.0/_r - 1.0);

  for (i=2;i<_n;++i)
    _pa[i] = _pa[i-1] + (_pa[i-1] - _pa[i-2]) + a;

  _pa[_n]=1.0;
}


// given a set of point p, that can represent a polyline, return the position of
// these points in parametric space 0-1.
/////////////////////////////////////////////////////////////////////
void
Curve::getParamPositions( const vector<Point3D>& _p,
                          vector<double>& _s )
{
 int i;
 double d;
 Point3D vec;
 int n = (int) _p.size();
 _s.resize (n);

 d = _s[0] = 0.0;

 for (i = 1; i < n; ++i)
 {
   double dist = _p[i].distance(_p[i-1]);
   d += _s[i] = dist;
 }
 for (i = 1; i < n; ++i)
   _s[i] = _s[i] / d + _s[i-1];
}


// given a line represented by LineStart-LineEnd points, return the distance d from this
// line to point Point. Obs: it returns 0 if the point p does not lie on the line
/////////////////////////////////////////////////////////////////////
int
Curve::distancePoint2Line( Point3D _LineStart,
                           Point3D _LineEnd,
                           Point3D _Point, 
                           Point3D& _Intersection,
                           double& _d )
{
  double LineMag2 = _LineStart.distance2(_LineEnd);

  if (LineMag2 == 0.0)
    return 0;

  Point3D v1 = _Point   - _LineStart;
  Point3D v2 = _LineEnd - _LineStart;

  double U = (v1.dotProd (v2)) / LineMag2;

  if( U < 0.0 || U > 1.0 )
    return 0; // closest point does not fall within the line segment

  _Intersection = _LineStart + U * v2;

  _d = _Point.distance(_Intersection);

  return 1;
}

// given a line segment represented by LineStart-LineEnd points, return
// the distance d from this line segment to the point Point.
// Obs: it returns 1 if t is inside the interval [0,1] and 0 otherwise.
/////////////////////////////////////////////////////////////////////
int
Curve::distancePoint2Segment( Point3D _LineStart,
                              Point3D _LineEnd,
                              Point3D _Point, 
                              Point3D& _Intersection,
                              double& _d )
{
  double LineMag2 = _LineStart.distance2(_LineEnd);

  if (LineMag2 == 0.0)
    return 0;

  Point3D v1 = _Point   - _LineStart;
  Point3D v2 = _LineEnd - _LineStart;

  double U = (v1.dotProd (v2)) / LineMag2;

  if( U < 0.0 || U > 1.0 ) // closest point does not fall within the line segment
  {
    // then the distance is the minimum of the distance to the bound points
    double d1 = _Point.distance( _LineStart );
    double d2 = _Point.distance( _LineEnd );
    if( d1 > d2 )
    {
      _Intersection = _LineEnd;
      _d = d2;
    }
    else
    {
      _Intersection = _LineStart;
      _d = d1;
    }
    return 0; 
  }

  _Intersection = _LineStart + U * v2;

  _d = _Point.distance( _Intersection );

  return 1;
}


//
/////////////////////////////////////////////////////////////////////
void 
Curve::getAdaptSubdvision( const double _maxAngle, 
                           vector<double>& _pa, bool adapt_only /*= false */)
{
  double dist;
  Point3D pt0, tg0, pt1, tg1;

  double len = getLength();
  set <double> param;

  // init begin and end parametric values
  this->get3dPoint (0.0, pt0, tg0);
  this->get3dPoint (1.0, pt1, tg1);
  dist = pt0.distance(pt1);
  param.insert(0.0);
  param.insert(1.0);

  // check adaptive refinement
  if (!checkAdaptiveEquation (_maxAngle, len, dist))
    adaptiveRefinement (_maxAngle, 0.0, 1.0, len, param);

  // fill vector
  set <double>::iterator itr;
  _pa.clear();
  // cout << "Refinamento adaptativo" << endl;
  for (itr = param.begin(); itr != param.end(); itr++)
  {
    double val = *itr;
    // fprintf (stdout, "%12.10lf ", val);
    _pa.push_back(val);
  }

  // consider max size of subdvision
  if (adapt_only == false)
  {
    vector <double> par_subdv;
    vector <Point3D> pts_subdv;
    double local_len, min_len = -1.0;

    // get points and parametric subdivision
    getSubdvision (pts_subdv);
    getParamSubdivision (par_subdv);

    // get min size
    int n = (int) par_subdv.size();
    int i, k;
    for (i = 0; i < n-1; ++i)
    {
      local_len = (par_subdv[i+1] - par_subdv[i]) * len;
      if (local_len < min_len || min_len == -1.0)
        min_len = local_len;
    }

    // compute new parametric subdvisions
    n = (int) _pa.size();
    for (i = 0; i < n-1; ++i)
    {
      double dt = _pa[i+1] - _pa[i];
      local_len = dt * len;
      if (local_len > min_len)
      {
        int local_subd = int (local_len / min_len);
        if (local_subd == 1)
          continue;
        
        for (k = 0; k < local_subd-1; k++)
          _pa.push_back (_pa[i] + dt / local_subd * (k+1));
      }
    }


    // reorder parametric subdvision
    sort (_pa.begin(), _pa.end());    
  }

  //int n = (int) _pa.size();
  //for (int i = 0; i < n; ++i)
  //{
  //  cout << _pa[i] << " " ;
  //}
  //cout << endl;
}

//
/////////////////////////////////////////////////////////////////////
bool 
Curve::checkAdaptiveEquation( double _maxAngle, double _li, double _di )
{
  if (_di == 0)
    return true;

  double ang_rad = _maxAngle * 3.1416 / 180.0;

  if ( (_li / _di) <= (ang_rad / (2 * sin (ang_rad * 0.5))) )
    return true;
  
  return false;
}

//
/////////////////////////////////////////////////////////////////////
void 
Curve::adaptiveRefinement( const double _maxAngle,
                           const double _par_i, 
                           const double _par_j, 
                           const double _len, 
                           set <double> &_param )
{
  double par_m = (_par_j+_par_i) * 0.5;
  Point3D pti, tgi, ptj, tgj, ptm, tgm;

  //    +-----------+------------------+--------------------+-----------+
  //    0.0         i                 mid                   j          1.0

  // add mid point
  _param.insert(par_m);

  // cout << "AdaptRef " << _par_i << " " << _par_j << " " << _len << endl;

  // eval points
  this->get3dPoint (_par_i, pti, tgi);
  this->get3dPoint (_par_j, ptj, tgj);
  this->get3dPoint ( par_m, ptm, tgm);

  // eval distances
  double dist_i_m, dist_m_j;
  dist_i_m = pti.distance(ptm);
  dist_m_j = ptm.distance(ptj);

  // check adaptive refinement for first segment
  if (!checkAdaptiveEquation (_maxAngle, _len*0.5, dist_i_m))
    adaptiveRefinement (_maxAngle, _par_i, par_m, _len*0.5, _param);

  // check adaptive refinement for second segment
  if (!checkAdaptiveEquation (_maxAngle, _len*0.5, dist_m_j))
    adaptiveRefinement (_maxAngle, par_m, _par_j, _len*0.5, _param);
}


//
/////////////////////////////////////////////////////////////////////
Curve*
Curve::join( Curve* _otherCurve,
             double _tol )
{
  //---------------------------------------------------
  // IMPLEMENTATIONs for the FUTURE: it should return a JOINED curve
  //---------------------------------------------------
  if( ((m_type==Curve::LINE)&&(_otherCurve->m_type==Curve::POLYLINE)) ||
      ((m_type==Curve::POLYLINE)&&(_otherCurve->m_type==Curve::LINE)) )
  {
    Line* line = NULL;
    PolyLine* poly = NULL;
    // defines which curve is the LINE and the POLYLINE
    if( m_type == Curve::LINE )
    {
      line = (Line *)this;
      poly = (PolyLine *)_otherCurve;
    }
    else
    {
      line = (Line *)_otherCurve;
      poly = (PolyLine *)this;
    }
    vector<Point3D> linePts = line->getCurvePoints();
    vector<Point3D> polyPts = poly->getCurvePoints();
    int npts = (int)polyPts.size();
    // starts the procedure to create the merged curve
    double tol = 1.0e-5;
    PolyLine* newCurve = new PolyLine;
    // tests if an endpoint has the same coordinates of a start point of the other
    int i;
    if( linePts[0].isAtSamePosition( polyPts[0], tol ) )
    {
      newCurve->addPoint( linePts[1] );
      for(i = 0; i < npts; ++i)
        newCurve->addPoint( polyPts[i] );
    }
    else if( linePts[0].isAtSamePosition( polyPts[npts-1], tol ) )
    {
      newCurve->addPoint( linePts[1] );
      for(i = npts-1; i >= 0; --i)
        newCurve->addPoint( polyPts[i] );
    }
    else if( linePts[1].isAtSamePosition( polyPts[0], tol ) )
    {
      newCurve->addPoint( linePts[0] );
      newCurve->addPoint( linePts[1] );
      for(i = 1; i < npts; ++i)
        newCurve->addPoint( polyPts[i] );
    }
    else if( linePts[1].isAtSamePosition( polyPts[npts-1], tol ) )
    {
      newCurve->addPoint( linePts[0] );
      newCurve->addPoint( linePts[1] );
      for(i = npts-2; i >= 0; --i)
        newCurve->addPoint( polyPts[i] );
    }
    newCurve->endAddPoint();
    return newCurve;
  }
  return 0;
}

